Non-Intestinal Microbial Signatures within Stool as Predictors of Cancer Immunotherapy Outcome

Introduction

Despite the increasing use of immunotherapy in treating advanced malignancies, clinical outcomes remain heterogeneous, with many patients exhibiting intrinsic resistance or eventual disease recurrence following immune checkpoint inhibitor (ICI) therapy [Robert et al., 2015]. This highlights the critical need for predictive biomarkers to optimize treatment selection and personalize therapeutic strategies. Accumulating evidence indicates a strong correlation between gut microbial composition and responsiveness to cancer immunotherapy [Peng et al., 2020]. While advancements in targeted therapies and immunotherapy have improved survival outcomes [Switzer et al., 2022], ICIs targeting cytotoxic T-lymphocyte-associated antigen 4 (CTLA-4) and programmed cell death protein 1 (PD-1) demonstrate limited efficacy in approximately 50% of patients with metastatic melanoma [Larkin et al., 2015] and carry the risk of immune-related adverse events [Robert et al., 2019].

Consequently, elucidating the mechanistic relationship between the gut microbiome and antitumor immunity has become a priority in translational cancer research. Indeed, the gut microbiome is increasingly recognized as a critical modulator of anti-tumor immunity and, consequently, response to cancer immunotherapy. Initial investigations, utilizing both preclinical models [Sivan et al., 2015] and clinical cohorts, particularly those with melanoma [Routy et al., 2018], have demonstrated a compelling link between gut microbial composition and immunotherapy efficacy. This link has been further supported by fecal microbiota transplantation studies in both gnotobiotic mice [Gopalakrishnan et al., 2018] and human patients [Baruch et al., 2021]. While much of the early research focused on melanoma [Tsakmaklis et al., 2023], current efforts are shifting towards identifying conserved microbial biomarkers and mechanisms applicable across diverse cancers [Gunjur et al., 2024; Lin et al., 2025], predicated on the hypothesis that fundamental interactions between the gut microbiome and the immune system are broadly relevant to immunotherapy response.

Despite numerous published studies, the biological mechanisms by which gut microbes regulate the immune system in the context of melanoma immunotherapy remain incompletely understood. A recent meta-analysis identified robust stool metagenomic markers associated with successful melanoma immunotherapy [Olekhnovich et al., 2023]. Notably, no taxonomic predictors of negative response were identified. Expanding on our earlier investigations, we employed genome-resolved metagenomics to gain a deeper understanding of the gut microbiota’s role in modulating immunotherapy response in malignant tumors. This analysis revealed a link between systemic metabolic depletion in the gut microbiome and decreased efficacy of immunotherapy in patients with metastatic melanoma [Zakharevich et al., 2024].

We expanded our application of genome-resolved metagenomics, analyzing 624 stool metagenomes from 11 published studies encompassing over five cancer types and representing diverse geographical regions across three continents. This analysis identified robust markers capable of differentiating patients based on immunotherapy response. These markers were validated using statistical analysis and machine learning techniques focused on logarithmic ratios of identified features. Notably, the analysis revealed non-intestinal traits such as oral- and food-derived bacteria within the marker list predominantly associated with negative treatment outcomes, potentially reflecting the previously observed depletion of gut microbial load associated with oral microbiota enrichment [Liao et al., 2024]. These findings underscore the importance of considering the metabolic landscape of the gut microbiome – and the surprising influence of non-intestinal microbial sources – when predicting and optimizing immunotherapy response, paving the way for the development of novel therapeutic strategies targeting microbiome-immune interactions.

References:

Materials and methods

Data acquisition and preprocessing

All 624 stool metagenomic samples from 11 published studies were downloaded from the NCBI/EBI databases using Kingfisher v0.4.1 [Woodcroft, 2022]. Data quality was assessed with FastQC v0.12.1 [Andrews et al., 2023], and raw reads were processed with fastp v0.23.4 to remove low-quality sequences [Chen et al., 2018]. Human-derived sequences were filtered out using HISAT2 v2.2 [Kim et al., 2019] with the GRCh37 human genome (Release 47) as a reference [Schneider et al, 2017].

Assembly and binning

Metagenomic assembly was performed using MEGAHIT v1.2.9 [Li et al., 2015], retaining contigs longer than 1000 bp. Filtered contigs were aligned to metagenomic reads using HISAT2 v2.2 [Kim et al., 2019]. Binning was conducted in two stages: (1) initial binning with MetaBAT2 v2.12.1 [Kang et al., 2019], MaxBin2 v2.2.7 [Wu et al., 2016], and Semibin2 v2.1.0 [Pan et al., 2023], followed by (2) refinement with DAS Tool v1.1.7 to improve bin quality [Sieber et al., 2018]. Bins were quality-checked with CheckM v1.2.3 [Parks et al., 2015] and dereplicated at 98% nucleotide identity using dRep v3.4.5 [Olm et al., 2017] to generate operational genomic units (OGUs) by analogy with operational taxonomic units (OTU). Zhu et al. (2022) first introduced this term, but our approach to representing OGEs offers a key advantage: OGEs are generated directly from the analyzed data, rather than through a mapping process. This more faithfully captures the core concept of the “operational unit”.

Taxonomic profiling and analysis

Taxonomic annotation of bins was performed using GTDB-Tk v2.1.1 [Chaumeil et al., 2022] with the GTDB r207 database [Parks et al., 2022]. Additional data processing utilized Samtools v1.17 [Li et al., 2009], BEDTools v2.31.0 [Quinlan et al., 2010], and BBMap v39.06 [Bushnell et al., 2014]. Finally, OGU abundance profiles were generated by mapping reads to the OGU catalog using HISAT2 v2.2, followed by processing with InStrain v1.9.0 [Olm et al., 2021].

Marker OGU discovery

The identification of OGUs associated with immunotherapy outcomes followed an established analytical framework from our previous works [Olekhnovich et al., 2023; Zakharevich et al., 2024]. Differential rankings were first performed using Songbird [Morton et al., 2019] to identify OGUs showing relative abundance variations between experimental groups, applying a conservative absolute differential value threshold of > 0.3. For candidate OGUs meeting this criterion, we subsequently calculated log-ratio abundances using Qurro [Fedarko et al., 2020] and determined statistical significance through Wilcoxon rank-sum tests implemented in the R statistical environment. The biomarker selection process incorporated stringent cross-validation criteria to ensure robust identification of microbial signatures. OGUs demonstrating consistent positive associations with therapeutic response across multiple datasets were retained as potential beneficial biomarkers, while any evidence of negative association with treatment outcome in any dataset resulted in automatic exclusion regardless of other positive associations.This approach enabled simultaneous identification of two clinically meaningful biomarker categories: microbial taxa positively correlated with successful immunotherapy outcomes and those associated with adverse therapeutic responses. The methodology emphasizes reproducibility through multi-dataset validation and maintains rigorous standards for biomarker qualification by requiring consistent directional effects across independent cohorts.

Markers testing

Determined R/NR biomarkers were used to calculate log ratios, followed by statistical assessment using the lmer function from the lmerTest package [Kuznetsova et al., 2017] in the R statistical environment. Post-hoc testing was performed using Tukey’s method with the multcomp package [Bretz et al., 2011], applying Bonferroni correction for multiple comparisons. Model residuals were tested for normality using the Shapiro-Wilk test (shapiro.test function in R). Additional testing of marker OGU included investigating the ability of calculated log ratios to separate experimental groups using logistic regression using the method of “leave-one group out” cross validation. Rahman et al. (2023) demonstrated the effectiveness of using log relations in machine learning models.Gene-set enrichment analysis (GSEA) implemented in GSEA function from clusterProfiler package [Xu et al., 2024] using to investigate statistical relationships of OGU marker sets across taxonomy and other attributes included food, body site or pathogens.GSEA analysis results were visualized using GseaVis package [Zhang et al., 2025].

Marker clustering with non-gut databases

The presence of non-gut microbes in the stool microbiome, such as those frequently detected on the oral means, has been associated with a degradation of the gut microbiota [Liao et al., 2024]. This phenomenon frequently results in adverse outcomes, including various illnesses. The identification of non-gut microbes was conducted using published microbial MAGs from body sites (oral, skin, vagina) [Pasolli et al., 2019] and food [Carlino et al., 2024]. The body site MAGs were retrieved from here. A total of 9,412 genomes were selected based on specific criteria, including those from the “Body Site” category, such as the airways, nasal cavity, oral cavity, skin, and vagina. Genomes classified as “Age Category” as either “newborn” or “child” were excluded from the analysis. The selected genomes were then clustered with our reconstructed catalog via dRep [Olm et al., 2017] with 95% nucleotide identity. The genome that exhibited a high degree of similarity with body MAG marks from the same source. Consequently, 293 MAGs were designated as body microbes. The food MAGs were retrieved from here, and the archive cFMD_mags.tar.gz contains the folder cFMD_mags_prok, which contains prokaryotic genomes. These 10,112 metagenome-assembled genomes (MAGs) were subsequently clustered via dRep with 95% nucleotide identity, yielding a total of 983 genomes. The subsequent source identification process was consistent with that employed for body MAGs. Of these, 91 were classified as food microbes.

Functional profiling of marker OGU

Functional profiles of OGUs were generated using MetaCerberus [Figueroa et al., 2024] with annotations from the Carbohydrate-Active EnZymes (CAZy) database [Drula et al., 2022] and the Kyoto Encyclopedia of Genes and Genomes (KEGG) [Kanehisa et al., 2025]. Statistical associations between functional features and experimental groups were assessed using logistic regression. Gene set enrichment analysis (GSEA) was applied to evaluate links between gene sets and immunotherapy response, with results visualized using the GseaVis package.

Additional analysis and visualization

For visualization and statistical analysis, we utilized the following R packages: ggplot2 [Wickham et al., 2016] for graphics creation, ggpubr [Kassambara, 2023] for publication-ready plots, hrbrthemes [Rudis, 2020] for enhanced themes, and RColorBrewer [Neuwirth, 2022] for color palette management. String data manipulation was performed using the stringr package [Wickham, 2022]. The data analysis report was generated with Quarto [Allaire et al., 2022], incorporating dynamic tables via DT [Xie et al., 2023] and enhanced download functionality using downloadthis [Sidi, 2023]. Genes were ranked for GSEA based on p-value multiplied by the sign of the odds ratio, where a negative sign indicates association with the NR group and a positive sign indicates association with the R group.

References

Results

Data overview

The analysis incorporated 624 stool metagenomes obtained from 11 independent external datasets. Patients were stratified by immunotherapy response into two groups: responders (R group, n = 362; 58.2%) and non-responders (NR group, n = 262; 41.8%). Response assessment followed RECIST 1.1 criteria, with the R group including patients showing complete response (CR), partial response (PR), or stable disease (SD) at 6-month follow-up, while the NR group comprised exclusively progressive disease (PD) cases. The study cohort received various immunotherapy regimens, including anti-PD1, anti-CTLA4, or combination therapies. Cancer type distribution revealed melanoma predominance (n = 456; 72.7%) [Frankel et al. 2017; Gopalakrishnan et al. 2018; Matson et al. 2018; Spencer et al. 2021; Lee et al. 2022; McCulloch et al. 2022; Tsakmaklis et al. 2023], followed by gastrointestinal (GI) cancers (n = 82; 13.1%) [Peng et al. 2020; Heshiki et al. 2020; Gunjur et al. 2024], non-small cell lung cancer (n = 15; 2.4%) [Liu et al. 2022; Heshiki et al. 2020], breast cancer (n = 4; 0.6%) [Heshiki et al. 2020], ovarian cancer (n = 2; 0.3%) [Heshiki et al. 2020], and other malignancies (n = 68; 10.8%) [Gunjur et al. 2024]. All samples were collected prior to treatment initiation to evaluate baseline microbiota status.

Figure 1

Quality control

Following quality filtering with fastp, aggregated quality metrics were compiled using MultiQC to summarize the FastQC reports (MultiQC are attached to the report below). Before filtering, approximately ~15bn reads (mean ± SD: 24 ± 15 mln reads per sample) were retained for downstream analysis, with an average of 34 ± 18% of reads per sample removed during quality control (see Table 2 for detailed filtering statistics).

MultiQC report

Taxonomic composition

A total of 3649 OGUs belonging to 12 bacterial phyla were detected in 627 samples. The most numerous of these were Firmicutes (n = 2443), Actinobacteriota (n = 590), Bacteroidota (n = 394), Proteobacteria (n = 140), Desulfobacterota (n = 27) and others (n = 55). OGU abundance profiles of stool samples presented in Table 3.

Marker OGUs discovery

Analysis of marker OGUs using Songbird and Qurro identified 637 significantly differentiated features between R and NR experimental groups. Of these, 298 were enriched in R, while 339 were enriched in NR. The most prominent microbial species enriched in responders included Blautia wexlerae (n = 49), Faecalibacterium prausnitzii (n = 27), Fusicatenibacter saccharivorans (n = 15), Gemmiger qucibialis (n = 14), and Blautia faecis (n = 8). Conversely, NR exhibited enrichment of distinct microbial taxa, including F. prausnitzii (n = 17), Veillonella parvula (n = 7), and Dialister invisus (n = 7). Although F. prausnitzii was detected in both groups, its significantly higher abundance in responders (27 vs. 17 OGUs) suggests potential strain-level or functional variations contributing to treatment outcome. The robust association of Blautia species, particularly B. wexlerae, with positive response highlights its potential as a predictive biomarker. Notably, the enrichment of oral bacteria Veillonella parvula and Dialister invisus in NR suggests a potential link between oral microbiome composition and immunotherapy efficacy. Detailed descriptions of marker OGUs are provided in Table 4. Songbird coefficients for each OGU marker are presented in Figure 2.

Statistical analysis utilizing linear mixed-effects models, accounting for dataset variables, revealed significant differences between R and NR groups (p < 0.0001). The inclusion of the dataset variable as a fixed effect was statistically significant (p < 0.0001), justifying the employed modeling approach. Residuals were normally distributed (Shapiro-Wilk test, p < 0.0001), and the observed effect size was large (Cohen’s d = 0.91). Log ratio values from the metagenomic samples are detailed in Table 5, with their distribution visualized across experimental groups in Figure 3. A predictive classifier, based on log ratios and logistic regression with Synthetic Minority Oversampling Technique (SMOTE) balancing, demonstrated a receiver operating characteristic (ROC) area under the curve (AUC) of 0.72 ± 0.12 (Figure 4), with other performance metrics detailed in Table 6 and Figure 5.

Saving 8 x 3 in image

Figure 2. The barplot depicts Songbird differential abundance coefficients (MEAN ± SD) for marker Operational Genomic Units (OGUs), with the x-axis representing individual OGUs and the y-axis showing the magnitude and direction of association coefficients.

Saving 4 x 3 in image

Figure 3. The boxplot shows log ratio (LR) values obtained using marker OGU. LR is highest in the R group (median = 3.0) and lowest in NR group (median = 1.4).

lmer test results


     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lmer(formula = log_ratio ~ response + (1 | dataset), data = df.log_ratio)

Linear Hypotheses:
            Estimate Std. Error z value Pr(>|z|)    
R - NR == 0   1.7005     0.1508   11.27   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- bonferroni method)

Dataset variable inportance

ANOVA-like table for random-effects: Single term deletions

Model:
log_ratio ~ response + (1 | dataset)
              npar  logLik    AIC    LRT Df Pr(>Chisq)    
<none>           4 -1274.6 2557.1                         
(1 | dataset)    3 -1283.1 2572.3 17.119  1  3.511e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residuals of the model normality testing


    Shapiro-Wilk normality test

data:  residuals(lmer.log_ratio)
W = 0.98021, p-value = 1.817e-07

Effect size estimation

Cohen's d |         95% CI
--------------------------
-0.91     | [-1.08, -0.74]

- Estimated using pooled SD.

Saving 4 x 4 in image

Figure 4. Receiver operating characteristic (ROC) curve for log regression performance using “leave-one group out” (LOGO) cross-validation.

Saving 4.5 x 3.5 in image

Figure 5. Distribution of log regression prediction metrics across datasets obtained by LOGO cross-validation.

Marker OGUs characterization

GSEA: taxonomy

Gene set enrichment analysis (GSEA) identified patterns of change related to microbial taxonomy (phylum, genus) and origin, distinguishing between microbes from non-intestinal body sites (mouth, skin, vagina) and environment (food). OGUs with 98% nucleotide identity were grouped by taxonomy or other characteristics to facilitate enrichment analysis within specific groups. Results indicated that the Bacteroidetes phylum was upregulated in the R group and Proteobacteria in the NR group (abs. NES > 2, GSEA adj.p < 0.01). Obtained results presented in Figure 6 Table 7. Bacteroidetes OGU included in core enrichment species presented in Table 8.Notably, enrichment of Proteobacteria OGUs included potential opportunistic species such as Klebsiella quasipneumoniae, Klebsiella michiganensis, Enterobacter ludwigii, Enterobacter kobei, Citrobacter youngae, Citrobacter portucalensis, Citrobacter freundii, and Haemophilus sp001815355 (see Table 9). At the genus level, Blautia, Bacteroides, Fusicatenibacter, Gemmiger, and Faecalibacterium were enriched in the R group, while Veillonella, Prevotella, Dysosmobacter, and Acetatifactor were enriched in the NR group (see Figure 6 and Table 10; GSEA adj.p < 0.01).

Clustering 3,816 OGUs from this study with 9,412 MAGs derived from different body sites and 10,112 food-derived MAGs revealed that 91 OGUs exhibited complementary to food-derived microbes (see Table 11), and 293 were complementary to those from non-intestinal sources, at 95% nucleotide similarity (Table 12), which corresponds to the species level. Based on the GSEA results, NR-linked OGUs were enriched by food-derived (see Figure 7 and Table 13) and oral (see Figure 8 and Table 14) microbial flora. Core food-derived OGUs comprised Bifidobacterium animalis, C. freundii, E. kobei, Enterococcus faecalis, Enterococcus faecium, Escherichia coli, K. michiganensis, K. quasipneumoniae, Lactobacillus gasseri, Limosilactobacillus reuteri, and Megasphaera elsdenii, showing substantial overlap with core enrichment patterns observed in Proteobacteria (Table 15). Core enrichment patterns within the oral cavity included OGU representatives of Veillonellales and Enterobacterales, specifically Anaeroglobus micronuciformis, Escherichia coli, Veillonella parvula, and Haemophilus sp001815355 (Table 16).

Saving 6 x 5.5 in image

Figure 6. GSEA results at the phylum level.

Saving 6 x 5.5 in image

Figure 6. GSEA results at the genus level.

GSEA: origin

Saving 4.5 x 3.5 in image

Figure 7. GSEA results at food-related OGUs.

Saving 4.5 x 3.5 in image

Figure 8. GSEA results at bodysite-related OGUs.

Functional analysis of marker OGU

Functional analysis of marker OGUs was performed using the MetaCerberus pipeline and CAZy and KEGG databases. Logistic regression was used to examine statistical relationships between functional groups and R/NR OGU sets. GSEA was then employed to identify links between CAZy/KEGG functional gene groups and R/NR OGU sets, utilizing p-values generated by the logistic regression.

Glycoside hydrolase (GH) functional groups were linked to the R OGU set (NES = 1.4, adj. p = 0.003). Obtained results presented in Table 17 and Figure 9. The core enrichment set included 30 glycoside hydrolase families, such as GH39 (alpha-L-iduronidase), GH42 (beta-galactosidase), GH43 (beta-xylosidase), GH51 (endoglucanase), and GH73 (lysozyme) (see Table 18). The most common microbiota genera associated with these enriched GH families were Blautia (31 genomes with >65% of core enrichment GH families), Fusicatenibacter (n = 11), Bacteroides (n=9), and Parabacteroides (n=5) (see Figure 9).

Analysis of enriched KEGG pathways revealed a significant association with the R OGU and NR OGU groups (abs. NES > 1, adj. p < 0.01). Obtained results presented in Table 19 and Figure 10. Specifically, the R OGU group exhibited enrichment for eight KEGG pathways, including: ko00340 (Histidine metabolism), ko01210 (2-Oxocarboxylic acid metabolism), ko00290 (Valine, leucine and isoleucine biosynthesis), ko00220 (Arginine biosynthesis), ko01230 (Biosynthesis of amino acids), ko00920 (Sulfur metabolism), ko00040 (Pentose and glucuronate interconversions), and ko01240 (Biosynthesis of cofactors). The NR OGU group, in contrast, showed enrichment for two KEGG pathways: ko00540 (Lipopolysaccharide biosynthesis) and ko02040 (Flagellar assembly). Further analysis revealed distinct bacterial genera associated with these functional pathways (see Figure 12). Blautia and Fusicatenibacter were primarily linked to the pathways enriched in the R OGU marker group, while Veillonella, Dialister, and Enterobacter were associated with ko00540 (Lipopolysaccharide biosynthesis) in the NR group. Similarly, Acetifactor, Roseburia, and Ruminoclostridium were linked to ko02040 (Flagellar assembly) enrichment in the NR group. These findings suggest potential mechanistic links between specific microbial taxa and the observed functional differences between the R OGU and NR OGU groups.

GSEA: CAZy

Saving 4.5 x 3.5 in image

Figure 9. GSEA results with CAZy database.

Saving 6 x 4 in image

Figure 10. Marker OGU contains CAZy core enrichment glycoside hydrolase within the genome.

GSEA: KEGG

Saving 7.5 x 5 in image

Figure 11. GSEA results with KEGG database.

Saving 12.5 x 18.5 in image

Figure 12. Marker OGU contains KEGG pathways core enrichment gene groups within the genome.

References:

Discussion

Our analysis reveals a distinct microbial signature associated with positive response to immunotherapy, prominently featuring enrichment of Blautia species, particularly Blautia wexlerae and Blautia faecis. This finding reinforces our previous meta-analyses identifying Blautia as a significant biomarker for favorable outcomes in melanoma immunotherapy [Olekhnovich et al., 2023; Zakharevich et al, 2024], and aligns with the growing recognition of this genus’s importance for human health [Liu et al., 2021]. The beneficial effects of Blautia may stem from its capacity to enhance gut barrier function through stimulation of mucus production and SCFA (acetate and propionate) generation, activating the Ffar2 receptor [Holmberg et al., 2024]. Furthermore, Blautia’s unique metabolic capabilities, including acetate production via the Wood-Ljungdahl pathway and vitamin B12 synthesis [Trischler et al., 2022; Zakharevich et al, 2024], may contribute to increased carbon availability and support microbial community function, ultimately bolstering host immunity. Interestingly, reduced Blautia abundance and acetate levels have been linked to diminished immune responses and accelerated cancer progression, particularly in contexts of chronic stress and depression [Ye et al., 2023]. These findings suggest that strategies to restore Blautia populations and acetate production hold promise as novel therapeutic interventions for enhancing immunotherapy efficacy and improving outcomes in patients.

Beyond Blautia, several other taxa – including species within the Bacteroides genus, Fusicatenibacter saccharivorans, and Gemmiger qucibialis – showed a positive correlation with response to immunotherapy, reinforcing our earlier findings identifying these species as biomarkers for favorable outcomes [Zakharevich et al, 2024]. While their specific contributions require further investigation, their co-occurrence with Blautia suggests a cooperative metabolic network centered on complex carbohydrate processing and SCFA production. Our data supports this, linking these species to key metabolic functions such as glycoside hydrolase activity, amino acid and vitamin synthesis, and sulfur metabolism. Elevated glycoside hydrolase levels likely reflect dietary fiber intake, facilitating the breakdown of polysaccharides into fermentable monosaccharides that fuel SCFA production – creating a positive feedback loop for gut health and immune modulation [Spencer et al., 2021]. Notably, Blautia’s capacity to generate acetate directly from CO2 via the Wood-Ljungdahl pathway expands the metabolic flexibility of the gut microbiome, allowing carbon fixation independent of dietary inputs. This highlights the dynamic interplay between the microbiome, host metabolism, and environmental carbon sources. While the link between amino acid/vitamin synthesis and immune function is well-established, our findings also reveal a connection to sulfur metabolism. Specifically, a decrease in H2S-consuming bacteria, including Blautia and Gemmiger, contributes to H2S accumulation, a characteristic observed in Crohn’s disease [Marcelino et al., 2023].

In contrast to the microbial signature associated with positive response, our analysis revealed a distinct composition in non-responding patients. Notably, the enrichment of taxa typically found in the oral microbiome – Veillonella parvula and Dialister invisus – was consistently observed in stool samples from those failing to respond to immunotherapy. We attribute this finding to a decrease in total microbial load in the large intestine [Liao et al., 2024], which correlates well with a decrease in functional capacity. Given the compositional nature of metagenomic data [Gloor et al., 2017], observed changes represent relative, not absolute, shifts in microbial abundance. These relative increases or decreases can arise from various scenarios, including an overall decline in the gut microbiota accompanied by a disproportionate loss in specific taxa, or an overall increase where one taxon expands more readily [Morton et al., 2019]. We hypothesize that a decline in the gut microbiota, coupled with the constant transit of oral microbes through the digestive system, creates a scenario where oral taxa serve as a relatively stable “reference” against which gut microbial representation is normalized. While we acknowledge the possibility that oral microbes may also actively colonize vacated niches. Interestingly, alongside increased oral taxa, we observed a cohort of bacteria frequently sourced from food, likely facilitated by a decrease in overall gut microbial load and the consistent supply of these bacteria through diet. While this included beneficial species commonly found in fermented foods and dairy products – such as Bifidobacterium animalis, Lactobacillus gasseri, and Limosilactobacillus reuteri – it also encompassed opportunistic pathogens like Citrobacter freundii, Enterobacter kobei, Enterococcus faecalis, Klebsiella michiganensis, and Klebsiella quasipneumoniae. This dual presence highlights the complex interplay between diet and microbiome composition. While fermented foods and probiotics can be beneficial, our findings underscore the critical importance of food source control. The introduction of opportunistic pathogens through ingested foods may contribute to dysbiosis and reduced immunotherapy efficacy. Therefore, both dietary composition and sourcing deserve careful consideration. This reveals a “double-edged sword” effect of dietary interventions: promoting a healthy microbiome through fiber intake is crucial, but equally important is minimizing the introduction of potentially pathogenic bacteria. Future studies should investigate the mechanisms underlying oral and food-borne microbial translocation and their impact on immunotherapy response. Targeted interventions focused on restoring gut barrier function, reducing microbial load, and ensuring pathogen-free food sources may offer promising strategies for enhancing immunotherapy efficacy in non-responding patients.

Interestingly, these observations align with the Anna Karenina principle applied to metagenomics [Zaneveld et al., 2017]. Healthy gut microbiomes are characterized by a conserved set of evolutionarily adapted bacteria – exhibiting a characteristic “sameness” (“all happy families are happy in the same way”) – whereas disruption of this core microbiome allows for increased detection of transient microbes from the oral cavity or environment, which constantly traverse the gut during ingestion. Consequently, greater degradation of the core gut microbiota results in a more variable microbial profile - “every unhappy family is unhappy in its own way”.

The contrasting microbial profiles observed in responders and non-responders highlight the complexity of host-microbiome interactions in modulating immunotherapy efficacy. While V. parvula and related oral commensals were significantly enriched in non-responders, the role of Faecalibacterium prausnitzii – a species often reported as a robust biomarker of response in other studies – proved more nuanced in our analysis. F. prausnitzii is frequently cited as a beneficial bacterium, known for its anti-inflammatory properties and its ability to produce butyrate, a short-chain fatty acid with enhanced immunotherapy effect [Gao et al., 2023; Bredon et 2024; Bredon et al., 2025]. However, our data, consistent with findings from our previous research [Olekhnovich et al., 2023; Zakharevich et al, 2024] and other researchers [Gunjur et al., 2024], reveal a more ambiguous association between F. prausnitzii abundance and immunotherapy response. We observed a non-consistent pattern, with F. prausnitzii levels demonstrating a variable behavior across both responder and non-responder groups. This seemingly contradictory finding could be attributable to several factors. Firstly, F. prausnitzii included multiple strains with potentially divergent functional roles. Different strains may exhibit varying levels of butyrate production or elicit distinct immune responses. Secondly, the abundance of F. prausnitzii may not be a reliable indicator of its activity. Functional studies assessing the metabolic capacity of F. prausnitzii within the gut environment are necessary to determine whether its presence correlates with genuine butyrate production and anti-inflammatory effects. Furthermore, the context of the broader microbial community likely influences the function of F. prausnitzii. In a dysbiotic gut, even the presence of F. prausnitzii may be insufficient to overcome the negative impact of other, pro-inflammatory species. The observed variable behavior may suggest that F. prausnitzii functions as a conditional biomarker – its benefit is only realized within a specific microbial context, or in the absence of other, competing factors. This finding highlights the limitations of not relying solely only on taxonomic abundance data for microbiome research. While identification of key bacterial taxa is a valuable starting point, a more complete understanding of microbial functions and complex interactions in the gut ecosystem is critical for the development of effective microbiome-directed therapies. Future studies should employ metagenomic and metabolomic approaches to characterize the functional capacity of F. prausnitzii in both responders and non-responders, and to identify the specific microbial interactions that determine its impact on immunotherapy response. Investigating the strain-level diversity of F. prausnitzii will be crucial for elucidating its complex role in modulating anti-tumor immunity.

References:

Additional information

Data availability

In this study, we used open access data from the NCBI/EBI Sequence Read Archives, identified by the following BioProjects accession numbers: PRJNA397906, PRJEB22893, PRJNA399742, PRJNA770295, PRJEB43119, PRJNA762360, PRJNA1011235, PRJNA615114, PRJNA866654, PRJNA494824, PRJEB49516. MAGs catalog assembly pipeline was described in https://github.com/JeniaOle13/Cancer_MAGs. All MAGs sequences were deposited in NCBI under accession PRJNA1196825. Songbird Qurro files for the QIIME2 viewer were available at Zenodo. Source code and Quarto report available at https://github.com/JeniaOle13/cancer-biomarkers.

Fundings

Financial support for this study was kindly provided by the Russian Science Foundation under Grant No. 22-75-10029, which can be accessed at https://rscf.ru/project/22-75-10029/.